#include "refiner.h"

/*** Refiner
 * constructor
 * ***/
#undef __FUNCT__
#define __FUNCT__ "Refiner"
Refiner::Refiner(MPI_Comm local_comm_, int local_color_, int num_refinements)
: local_comm(local_comm_), local_color(local_color_), n_sims(num_refinements)
{
  solvers = new Solver*[n_sims];
}

/*** ~Refiner
* destructor
* ***/
#undef __FUNCT__
#define __FUNCT__ "~Refiner"
Refiner::~Refiner()
{
  for(int k = 0; k < n_sims; k++) {
    delete solvers[k];
  }
  delete[] solvers;
}

/*** run
 * run the simulation
 * ***/
#undef __FUNCT__
#define __FUNCT__ "run"
PetscErrorCode Refiner::run()
{
  Options               *opt = new Options();
  opt->read_options();
  
  for(int k = 0; k < n_sims; k++) {
#ifdef VERBOSE
    PetscPrintf(local_comm, "[INFO] ********** GROUP %d: Starting simulation %d of %d **********\n", local_color, k, n_sims);
#endif
    solvers[k] = new Solver(opt, local_comm, local_color);
    ierr = solvers[k]->solve(); CHKERRQ(ierr);
    
    opt->Nx *= 2;
    opt->Ny *= 2;
  }
  
  delete opt;
  
  return 0;
}

/*** save_error
 * save errors data and plot them
 * ***/
#undef __FUNCT__
#define __FUNCT__ "save_error"
PetscErrorCode Refiner::save_error( PlotType plot )
{
  if( solvers[0]->loc_rank != 0 ) return 0;
  
  FILE                  *f;
  char                  fname[PETSC_MAX_PATH_LEN], dat_fname[PETSC_MAX_PATH_LEN], pdf_fname[PETSC_MAX_PATH_LEN], command[PETSC_MAX_PATH_LEN];
  
  sprintf(fname, "%s%s_error_%.2f_nproc_%d_sim_%d", PATH_TO_DATA, solvers[0]->name, solvers[0]->T, solvers[0]->comm_size, local_color);
  strcpy(dat_fname, fname); strcat(dat_fname, ".dat");
  strcpy(pdf_fname, fname); strcat(pdf_fname, ".pdf");
  
#if VERBOSE
  PetscPrintf(PETSC_COMM_SELF, "[INFO] (Process %d of %d, group %d) Exporting error data to \"%s\"...\n",
              solvers[0]->loc_rank, solvers[0]->comm_size, local_color, dat_fname);
#endif
  
  f = fopen(dat_fname, "w");
  for ( int k = 0; k < n_sims - 1; k++ )
  {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%G %G %G %G\n", solvers[k]->dx, l2.at(k), l1.at(k), linf.at(k)); CHKERRQ(ierr);
  }
  fclose(f);
  
  if( plot != NONE ) {
    integer usepdf;
    if( plot == PDF) usepdf = 1;
    else usepdf = 0;
    sprintf(command, "gnuplot -e \"inname='%s';outname='%s';usepdf=%d\" './tools/errors.gnuplot'",
            dat_fname, pdf_fname, usepdf);
#if VERBOSE
    PetscPrintf(PETSC_COMM_SELF, "[INFO] (Process %d of %d) Running \"%s\"...\n",
                solvers[0]->loc_rank, solvers[0]->comm_size, command);
#endif
    system(command);
  }
  
  return 0;
}

/*** save_runtime
 * save runtime data and plot them
 * ***/
#undef __FUNCT__
#define __FUNCT__ "save_runtime"
PetscErrorCode Refiner::save_runtime( PlotType plot )
{
  if( solvers[0]->loc_rank != 0 ) return 0;
  
  FILE                  *f;
  char                  fname[PETSC_MAX_PATH_LEN], dat_fname[PETSC_MAX_PATH_LEN], pdf_fname[PETSC_MAX_PATH_LEN], command[PETSC_MAX_PATH_LEN];
  
  sprintf(fname, "%s%s_runtime_%.2f_nproc_%d_sim_%d", PATH_TO_DATA, solvers[0]->name, solvers[0]->T, solvers[0]->comm_size, local_color);
  strcpy(dat_fname, fname); strcat(dat_fname, ".dat");
  strcpy(pdf_fname, fname); strcat(pdf_fname, ".pdf");
  
#if VERBOSE
  PetscPrintf(PETSC_COMM_SELF, "[INFO] (Process %d of %d) Exporting runtime data to \"%s\"...\n",
              solvers[0]->loc_rank, solvers[0]->comm_size, dat_fname);
#endif
  
  f = fopen(dat_fname, "w");
  for ( int k = 0; k < n_sims; k++ )
  {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%G %G\n", solvers[k]->dx, solvers[k]->elapsed); CHKERRQ(ierr);
  }
  fclose(f);
  
  if( plot != NONE ) {
    integer usepdf;
    if( plot == PDF) usepdf = 1;
    else usepdf = 0;
    sprintf(command, "gnuplot -e \"inname='%s';outname='%s';usepdf=%d\" './tools/runtime.gnuplot'",
            dat_fname, pdf_fname, usepdf);
#if VERBOSE
    PetscPrintf(PETSC_COMM_SELF, "[INFO] (Process %d of %d) Running \"%s\"...\n",
                solvers[0]->loc_rank, solvers[0]->comm_size, command);
#endif
    system(command);
  }
  
  return 0;
}

/*** save_report
 * save report of simulation
 * ***/
#undef __FUNCT__
#define __FUNCT__ "save_report"
PetscErrorCode Refiner::save_report()
{
  if( solvers[0]->loc_rank != 0 ) return 0;
  
  FILE                  *f;
  char                  fname[PETSC_MAX_PATH_LEN], txt_fname[PETSC_MAX_PATH_LEN];
  
  real                  rate, E1, E2;
  
  sprintf(fname, "%s%s_report_%.2f_nproc_%d_sim_%d", PATH_TO_DATA, solvers[0]->name, solvers[0]->T, solvers[0]->comm_size, local_color);
  
  strcpy(txt_fname, fname); strcat(txt_fname, ".log");
  
  #if VERBOSE
  PetscPrintf(PETSC_COMM_SELF, "[INFO] (Process %d of %d, group %d) Exporting report to \"%s\"...\n",
              solvers[0]->loc_rank, solvers[0]->comm_size, local_color, txt_fname);
  #endif
  
  f = fopen(txt_fname, "w");
  
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "*** Report generated the date %d by process %d of %d. ***\n\n", 0,
                      solvers[0]->loc_rank,solvers[0]->comm_size); CHKERRQ(ierr);
  
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "[1] GENERAL INFORMATION\n"); CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "     Name:              %s\n", solvers[0]->name); CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "     Time stepping:     %s\n", solvers[0]->adaptive_ts ==  true ? "adaptive" : "constant"); CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "     Lambda:            %.2f\n", solvers[0]->lambda); CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "     Final time:        %.2f\n", solvers[0]->T); CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "     Domain:            [%.4f, %.4f] x [%.4f, %.4f]\n",
                      solvers[0]->left, solvers[0]->right, solvers[0]->bottom, solvers[0]->top); CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "     grid sizes:        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%dx%d", solvers[k]->Nx, solvers[k]->Ny); CHKERRQ(ierr);
    if( k != n_sims - 1 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, ", "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n\n"); CHKERRQ(ierr); }
  }
  
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "[2] CONVERGENCE AND RUNTIME INFORMATION\n"); CHKERRQ(ierr);
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Mesh (nodes)                     "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%4d x%4d", solvers[k]->Nx, solvers[k]->Ny); CHKERRQ(ierr);
    if( k != n_sims - 1 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "  "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Mesh (dx)                        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", solvers[k]->dx); CHKERRQ(ierr);
    if( k != n_sims - 1 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Mesh (dy)                        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", solvers[k]->dy); CHKERRQ(ierr);
    if( k != n_sims - 1 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Runtime                          "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", solvers[k]->elapsed); CHKERRQ(ierr);
    if( k != n_sims - 1 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   [Rate]                        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 1; k++ ) {
    rate = log(solvers[k+1]->elapsed / solvers[k]->elapsed) / log(solvers[k+1]->dx / solvers[k]->dx);
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, " (%.2f) ", rate); CHKERRQ(ierr);
    if( k != n_sims - 2 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Error (L2)                       "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 1; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", l2.at(k)); CHKERRQ(ierr);
    if( k != n_sims - 2 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   [Rate]                        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 2; k++ ) {
    rate = log(l2.at(k+1) / l2.at(k)) / log(solvers[k+1]->dx / solvers[k]->dx);
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, " (%.2f) ", rate); CHKERRQ(ierr);
    if( k != n_sims - 3 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Error (L1)                       "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 1; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", l1.at(k)); CHKERRQ(ierr);
    if( k != n_sims - 2 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   [Rate]                        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 2; k++ ) {
    rate = log(l1.at(k+1) / l1.at(k)) / log(solvers[k+1]->dx / solvers[k]->dx);
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, " (%.2f) ", rate); CHKERRQ(ierr);
    if( k != n_sims - 3 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Error (Linf)                     "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 1; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", linf.at(k)); CHKERRQ(ierr);
    if( k != n_sims - 2 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   [Rate]                        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 2; k++ ) {
    rate = log(linf.at(k+1) / linf.at(k)) / log(solvers[k+1]->dx / solvers[k]->dx);
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, " (%.2f) ", rate); CHKERRQ(ierr);
    if( k != n_sims - 3 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Total Energy                     "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", solvers[k]->energy.back()); CHKERRQ(ierr);
    if( k != n_sims - 1 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   [Rate]                        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 1; k++ ) {
    E1 = solvers[k]->energy.back();
    E2 = solvers[k+1]->energy.back();
    rate = log(E2 / E1) / log(solvers[k+1]->dx / solvers[k]->dx);
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, " (%.2f) ", rate); CHKERRQ(ierr);
    if( k != n_sims - 2 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "Total Energy variation (dE)      "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", abs( solvers[k]->energy.back() - solvers[k]->energy.front() )); CHKERRQ(ierr);
    if( k != n_sims - 1 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   [Rate]                        "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims - 1; k++ ) {
    E1 = abs( solvers[k]->energy.back() - solvers[k]->energy.front() );
    E2 = abs( solvers[k+1]->energy.back() - solvers[k+1]->energy.front() );
    rate = log(E2 / E1) / log(solvers[k+1]->dx / solvers[k]->dx);
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, " (%.2f) ", rate); CHKERRQ(ierr);
    if( k != n_sims - 2 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  ierr = PetscFPrintf(PETSC_COMM_SELF, f, "sup                              "); CHKERRQ(ierr);
  for( int k = 0; k < n_sims; k++ ) {
    ierr = PetscFPrintf(PETSC_COMM_SELF, f, "%.2e ", solvers[k]->sup); CHKERRQ(ierr);
    if( k != n_sims - 1 ) { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "   "); CHKERRQ(ierr); }
    else { ierr = PetscFPrintf(PETSC_COMM_SELF, f, "\n"); CHKERRQ(ierr); }
  }
  // MiSSING TODO Linf divergence, sup_u and sup_v
  fclose(f);
  
  return 0;
}

/*** compute_error
 * parse each solution and compute error
 * ***/
#undef __FUNCT__
#define __FUNCT__ "compute_error"
PetscErrorCode Refiner::compute_error()
{
  if( solvers[0]->loc_rank != 0 ) return 0;
  
  scalar                *V1, *V2;
  integer               I_u, I_v, J_u_sw, J_u_se, J_u_nw, J_u_ne, J_v_sw, J_v_se, J_v_nw, J_v_ne;
  real                  err_u, err_v, sup_u, sup_v, l_1, l_2, temp;
  
#ifdef VERBOSE
  PetscPrintf(PETSC_COMM_SELF, "[INFO] Computing error...\n");
#endif
  
  for( int k = 0; k < n_sims - 1; k++ ) {
    sup_u = 0;
    sup_v = 0;
    l_1 = 0;
    l_2 = 0;
    
    if(solvers[k]->Nx * 2 != solvers[k+1]->Nx && solvers[k]->Ny * 2 != solvers[k+1]->Ny) {
      ierr = PetscPrintf(PETSC_COMM_SELF, "[ERROR] Invalid mesh refinement (%dx%d vs %dx%d)!\n",
                         solvers[k]->Nx, solvers[k]->Ny, solvers[k+1]->Nx, solvers[k+1]->Ny); CHKERRQ(ierr);
    }
  
    err.push_back( new scalar*[solvers[k]->Ny] );
    for( int j = 0; j < solvers[k]->Ny; j++ ) {
      err.at(k)[j] = new scalar[solvers[k]->Nx];
    }
  
    ierr = VecGetArray(solvers[k]->U_collect, &V1); CHKERRQ(ierr);
    ierr = VecGetArray(solvers[k+1]->U_collect, &V2); CHKERRQ(ierr);
    for ( int j = 0; j < solvers[k]->Ny; j++ )
    {
      for ( int i = 0; i < solvers[k]->Nx; i++ )
      {
        I_u = (i + j * solvers[k]->Nx) * 2;
        I_v = I_u + 1;
        J_u_sw = (i * 2 + j * solvers[k+1]->Nx * 2) * 2;
        J_u_se = J_u_sw + 2;
        J_u_nw = J_u_sw + 2 * solvers[k+1]->Nx;
        J_u_ne = J_u_nw + 2;
        J_v_sw = J_u_sw + 1;
        J_v_se = J_u_se + 1;
        J_v_nw = J_u_nw + 1;
        J_v_ne = J_u_ne + 1;
      
        err_u = V1[I_u] - 0.25 * ( V2[J_u_sw] + V2[J_u_se] + V2[J_u_nw] + V2[J_u_ne] );
        err_v = V1[I_v] - 0.25 * ( V2[J_v_sw] + V2[J_v_se] + V2[J_v_nw] + V2[J_v_ne] );
      
        err.at(k)[j][i] = err_u;
        
        //if( k == 1 ) PetscPrintf(PETSC_COMM_SELF, "%f ", err_u);
      
        temp = err_u * err_u + err_v * err_v;
      
        sup_u = max(abs(err_u),sup_u);
        sup_v = max(abs(err_v),sup_v);
      
        l_1 += sqrt(temp) * solvers[k]->dx * solvers[k]->dy / solvers[k]->area;
        l_2 += temp * solvers[k]->dx * solvers[k]->dy;
      }
      //if( k == 1 ) PetscPrintf(PETSC_COMM_SELF, "\n");
    }
    ierr = VecRestoreArray(solvers[k]->U_collect, &V1); CHKERRQ(ierr);
    ierr = VecRestoreArray(solvers[k+1]->U_collect, &V2); CHKERRQ(ierr);
  
    //PetscPrintf(PETSC_COMM_SELF, "sup_u = %f, sup_v = %f, l1 = %f, l2 = %f\n", sup_u, sup_v, l_1, l_2);

    l_2 = sqrt( l_2 ) / solvers[k]->area;
    linf.push_back( max(sup_u, sup_v) );
    l1.push_back( l_1 );
    l2.push_back( l_2 );
  }
  
  return 0;
}
